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ABSTRACT 

It has previously been shown that varying the numerical timestep during a 
symplectic orbital integration leads to a random walk in energy and angular 
momentum, destroying the phase space-conserving property of symplectic inte- 
grators. Here we show that when altering the timestep symplectic correctors 
can be used to reduce this error to a negligible level. Furthermore, these correc- 
tors can also be employed to avoid a large error introduction when changing the 
Hamiltonian's partitioning. We have constructed a numerical integrator using 
this technique that is nearly as accurate as widely used fixed-step routines. In 
addition, our algorithm is drastically faster for integrations of highly eccentricitic, 
large semimajor axis orbits, such as those found in the Oort Cloud. 

Subject headings: methods: numerical, Oort Cloud, solar system: general 
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1. Background 



In the past twenty years, mixed variable symplectic integrators have revolutionized 
dynamical modeling in planetary astronomy. These integrators work by splitting a system 's 



Hamiltonian into two s eparate Hamiltonians, Hx ep and Hj nt (I Wisdom fc Holman 



1991 



Saha fc Tremaine 



1992). Hxep is the Hamiltonian for a Keplerian orbit about a large 



central point mass, while Hi nt , or the interaction Hamiltonian, is the difference between 
the real Hamiltonian and H^ep- These two Hamiltonians are then integrated in a leapfrog 
algorithm, where a particle is first evolved under the equations of motion of Hx e p for a half 
time step (t/2), then under the equations of motion of Hj nt for a full time step (r), and 
finally under Hx ep again for t/2. It can be shown that thi s is a second-order approximation 



of th e evolution of the real Hamiltonian for a full step, r (ISaha fc Tremaine 



1992 



Yoshida 



19931 ). (It should be noted that the integration kernel in many popular simulation codes 
actually integrates H Int for half-steps and H Kep for a full step. This is also a second order 
symplectic algorithm.) 

Mixed variable symplectic methods have two major advantages over previous 
integration techniques. First, symplectic methods can enable much larger integration 
step sizes compared to other integration techniques if the system being modeled is nearly 
integrable (in the case of planetary systems, if Hxe P » Hj nt ). For most simulations of 
a planetary system, Hx e p is a fairly close approximation to the real Hamiltonian of the 
system. As a result, H Int can be viewed as a perturbation, e, to the Kepler problem, 
resulting in an integration error of order er 2 . Thus, the appropriate step size is simply 
a function of the perturbation magnitude. For typical solar system integ rations, roughly 



10-20 steps per orbit are suitable for an accurate symplectic integration (IMorbidelli 



Wisdom &: Holman 



2002 



19921 ). The second advantage of symplectic techniques is that the 



numerical errors for a system's integrals of motion can be expressed in Hamiltonian form 
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( jSaha &: Tremaindll992l ; lYoshidalll993l ). Hence, for sufficiently small step sizes the errors 
will not grow with time. 

During a symplectic integration, the numerical Hamiltonian conserves quantities 
analogous to integrals of motion, but the traditional integrals of motion (e.g. energy) of the 
numerical system actually oscillate about those of the real Hamiltonian. Furthermore, this 
oscillation has a period and ampli t ude determined by the chosen timestep and Hamiltonian 



partition ( ISaha fc Tremaine 



19921 ) . 



Wisdom et al 



( 119961 ) show that this oscillation takes 



the form of an infinite series and can be well-represented by its leading terms. Fortunately, 
these terms can be solved for and applied as a symplectic corrector at any point during 
a symplectic integration to reduce numerical error to order e 2 r 2 . Thus, when e is small 
Wisdom-Holman style symplectic methods have the potential to be much more accurate 
than their leapfrog form initially implies. 

As with all methods, symplectic integrators also suffer from a few limitations. To remain 
symplectic, several numerical parameters must remain fixed throughout an integration. One 
of these is the integration timestep. The numerical Hamiltonian is a function of the chosen 
timestep, and thus changing this step will generate a different numerical Hamiltonian. 
Unless the integrals of motion of the numerical and real Hamiltonian happen to agree 
exactly when the timestep is changed, these quantities will begin oscillating about values 
different from the real system's. An example of this type of error introduction is shown 
in the top panel of Figure [U, where we integrate a particle on an orbit with semi-major 
axis a = 300 AU and perihelion distance q = 100 AU for 30,000 yrs in the potential of 
the Sun and four giant planets. This particle never comes close to the planetary region. 
Consequently, the particle receives only incredibly small energy kicks from the planets 
during perihelion passages, and the orbital energy variations shown in the figure are almost 
purely due to numerical error. Initially, we integrate with a timestep of 0.548 yrs (200 days), 
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and even though the energy oscillations are large, they are also bounded. At t =10,000 yrs 
we decrease the timestep by half, and although the energy error decreases substantially 
with the smaller step size, it is clear that the numerical Hamiltonian is now modeling a 
system different from the original one. Finally, this same error is demonstrated again when 
we switch back to the original timestep at t =20,000 yrs and one can see that the numerical 
Hamiltonian drifts even further from the original system. 

Another numerical parameter that cannot be altered without compromising an 
algorithm's symplectic properties is the way the Hamiltonian is partitioned. The choice 
of canonical coordinates that are used in a simulation typically determines how the 
Hamiltonian is split between Hx e p and Hi nt . Changing this splitting will result in an 
error analogous to a time step change. Once again, we illustrate this effect in the bottom 
panel of Figure [TJ here we integrate the same particle as in the top panel, but in this 
case the timestep is held constant for the entire integration. Initially, the integrator uses 
a symplectic algorithm that performs Keplerian drifts about the Sun (in heliocentric 
canonical coordinates). Then, at t =10,000 yrs the Hamiltonian is reexpressed in canonical 
Jacobian coordinates (which are barycentric for massless test particles). Now, the drifts 
are performed about the solar system barycenter. Because Hj nt is so much smaller at large 
distances with this Hamiltonian formulation, the energy oscillations are greatly diminished. 
As in the top panel, however, we see that the numerical Hamiltonian is now tracking the 
behavior of a system that is different from our initial integration. Furthermore, switching 
back to the heliocentric algorithm at t =20,000 yrs shows that the particle is now on a 
slightly different orbit than it started with, and additional error is introduced. 

It should be noted that changing either the canonical coordinates or the timestep a 
handful of times during an integration will not typically invalidate a simulation. After 
a single change, the resulting numerical Hamiltonian and error magnitude will still be 
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Fig. 1. — Plots of numerical energy error vs. time for an orbit with a =300 AU, q =100 AU, 
and i =45°. Top Panel: A heliocentric symplectic integrator is used for this integration. 
The solid line corresponds to periods when the integration step was set to 0.548 yrs, and the 
dotted line marks periods using a 0.27-yr step. Bottom Panel: A symplectic integrator is 
used for this integration. The time step is held at 0.548 yrs for the entire integration. The 
solid line corresponds to periods when the Keplerian drift is performed about the Sun, and 
the dotted line corresponds to a Hamiltonian partition using a barycentric drift. 
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very clos e to that of the initial sim ulation. In the widely used SWIFT RMVS3 simulation 
package ( Levison Sz Duncan] 



19941 ). both the timestep and canonical coordinate system are 
changed during encounters between massless test particles and massive bodies. However, 
because these encounters are highly chaotic and occur relativel y infrequently, the small 



errors introduced do not alter the actual dynamics of particles (jLevison fc Duncan 



1994 



Migliorini et al. 



19971 ). However, when timestep and coordinate changes are made many 



times, each individual transition manifests itself as a step in a random walk away from the 
true integrals of motion of the system, and in certain situations this random walk can be 
the driving force behind the dynamics of orbits rather than the real Hamiltonian of the 
system. 

While the above limitations of symplectic routines are inconsequential for many 
scenarios, there are classes of problems when they impede accurate modeling. One such 
class is simulating planetary accretion. During close encounters between massive bodies in 
accretion simulations, Hj nt temporarily dominates Hx e p, and its treatment as a perturbation 
to Hxep is no longer valid, which causes the integrator to fail. Severa l appr oaches have 
been devised to circumvent this problem. For example, iDuncan et al.l ( 119981 ) decompose 
planetary forces into a series of forces acting over different ranges and sampled at different 
rates. Only during an encounter is it necessary to sample the high-frequency close-range 
forces sinc e they will be zero when massive bodies are far from one another. Alternatively, 
Chambers! (119991 ) temporarily switches to a Bulirsch-Stoer technique to integrate an Hx e p 
expression that absorbs encounter terms when they become too large. Although this second 
technique is not strictly symplectic, the Bulirsch-Stoer integratio n can be performed to 



machine precision resulting in no additional accuracy loss. Lastly, 



McNeil k Nelson! (J2009|) 



combin e these methods with a leapfrog approach assigning different fixed steps to different 



bodies ( iSaha fc Tremaine 



19941 ) to obtain a quasi-symplectic algorithm that successfully 



handles encounters and accurately adjusts the integration step for different radial zones. 
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Oort Cloud dynamics is another regime hampered by integrator limitations. When 
simulating the Oort Cloud, it is necessary to accurately integrate extremely eccentric 
test particle orbits. Consequently, a symplectic integrator whose H Kep is centered on the 
solar system baryce nter cannot be used sinc e this approximation fails during very close 



perihelion passages ( iRauch &: Holman 



19991 ). For this reason, an algorithm such as the 



SWIFT RM VS3 integrato r , whic h performs Keplerian drifts about the Sun, is typically used 



instead, e.g. 



Pones et al. 



( 12004 ) . However, this algorithm is quite inefficient for distant 
cometary orbits. Oort Cloud bodies orbit the Sun on Myr-timescales, but the simulation 
step size must be held below ~0.5 yrs (~ 5% of Jupiter's period) in order to accurately 
model the perturbations of the giant planets. This limitation even holds for orbits that do 
not approach the planetary region because integrations in the heliocentric frame always 
contain an indirect acceleration term due to planetary perturbations on the Sun, which 
does not fall off with particle distance. 

We present an accurate integration package called SCATR (Symplectically Corrected 
Adaptive Timestep Routine) that symplectically integrates test particles in a barycentric 
frame with large timesteps while they are far from the planetary region of the Solar System. 
When test bodies approach the planetary region, the integrator switches to heliocentric 
coordinates and uses a much smaller step size. To do this smoothly, the algorithm uses a 
symplectic corrector to drastically reduce the numerical error from time step changes as 
well as differences in the Hamiltonian partition. Although the timestep transitions and 
Hamiltonian reformulations in our algorithm are not symplectic, we show that the net 
numerical errors accrued by these transitions over the history of the Solar System are still 
comparable to, or even smaller than, the numerical error of an uncorrected symplectic 
integration. The benefit of this additional small numerical error is that the computing time 



of distant orbital integrations can be decreased by up to a fac tor of ~ 30. While not as 



elegant as the adaptive step scheme of 



McNeil fc Nelson! (120091 ). our algorithm is simple and 
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efficient. As such, it is well-suited for integrations focused purely on test particle dynamics 
such as Oort Cloud studies. 



2. The Algorithm 



Our numerical code, which is based on the routines found in SWIFT's RMVS3 



( ILevison fc Duncan 



19941 ). uses the following symplectic integration operator: 



S(t) =Z(t/2)o/C(t)oJ(t/2) (1 

where the operator X(t) corresponds to evolution under the equations of motion of Hj nt 
for a time t, K corresponds to evolution under the equations of motion of Hx ep , and £{r) 
corresponds to evolution of the entire numerical Hamiltonian over one step, r. Unlike 
SWIFT, ours is a hybrid c ode that only uses the re gularized mixed variable symplectic 



(RMVS) mapping scheme (ILevison &: Duncan 



19941 ) when test particles are close to the 



planetary region. When test particles are far from the planetary region our code uses the 
traditional Wisdom-Holman (WH) sym plectic mapping, which pe rforms Keplerian drifts 



about the barycenter for test particles (jWisdom fe Holman 



1991 



). These two mapping 



schemes differ by the assumed mass about which a particle is orbiting. For the WH scheme, 
it is assumed that particles are following nearly Keplerian orbits about the center-of-mass of 
the solar system. For test particles under the gravitational influence of the Sun and planets 



H 



Kep — 



GM 



(2) 



in this mapping scheme, where Vf, is the barycentric particle velocity, M is the solar system 
mass, and R is the particle distance to the barycenter. In this scheme the perturbing 
Hamiltonian is 

GM GM S ^GM p 
riint = — ^ > v — - — (a J 



R 



R.< 
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where Ms is the mass of the Sun, Rs is the distance from the particle to the Sun, the M p 
terms are the planetary masses, and the r p terms are the particle-planet distances for each 
of the planets in the system. 

Alternatively, the RMVS mapping scheme assumes that particles are following nearly 
Keplerian orbits about the Sun rather than the barycenter. Consequently, H^ep and Hj nt 
take on different forms in this map. For a test particle, the Keplerian Hamiltonian is now 



where vs is the heliocentric velocity of the test particle. For this same particle, the 
interaction Hamiltonian is now 



In this expression, the last summation represents the "indirect" acceleration experienced by 
the particle due to the acceleration of the Sun (the RMVS coordinate system origin) by the 
planets orbiting it. 

This map switching in our code is done for two reasons. First, unlike the WH map, 
RMVS mapping ensures an accurate integration through perihelion passage for any orbital 
eccentricity since the Keplerian drift is performed about the Sun. Second, while Hj nt for 
the RMVS method will always contain a large indirect acceleration term due to planetary 
perturbations on the Sun, the WH map approaches a pure Keplerian drift at large distances 
since Hj nt approaches zero there. If we can devise a method to smoothly switch timesteps 
and maps (which we detail below), computing efficiency may be greatly enhanced by using 
much larger integration steps for particles at large distances. 

As shown in Figure [TJ varying the Hamiltonian partition or integration timestep of 
a symplectic integrator introduces a large amount of error due to spurious oscillations of 
the system's integrals of motion under the numerical Hamiltonian. However, symplectic 




Rs 



(4) 
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19961 ) . and outfitting 



correctors remove a large degree of this oscillation (I Wisdom et al. 
our code with a corrector enables us to acc urately switch t imesteps and symplectic maps 
far from the Sun. Following the notation of I Wisdom! (120061 ). symplectic correctors can be 



expressed in terms of K and X. First, let 



X(ar, br) = K{ar) o X(6r) o /C(-ar). (6) 
where a and b are defined coefficients for a corrector of a given order. Next let 

Z(ar, br) = X(ar, br) o X(-ar, -br). (7) 

Now, an (n + l)-order corrector can be written as 

C(t) = Z(aiT, bir) o Z(a 2 r, b 2 r) o ... o Z(a n r, b n r). (8) 

In addition, an inverse corrector to convert from the real Hamiltonian to the numerical one 
can be expressed as 

C _1 (t) = Z(a n r, -b n r) o ... o Z(a 2 r, -b 2 r) o Z(a 1 r, -b x r). (9) 

Thus, a corrected symplectic algorithm that takes m steps will have the following form: 

E\mr) = C-\r) o [X(r/2) o JC(t) o X(r/2)] m o C(r). (10) 

or 

S'(mr) =C- 1 (T)o[S(r)] m oC(T) (11) 

Note that the corrector and its inverse only need to be applied at the very beginning and 
at the end of the sequence of m steps since they will cancel each other out in between 
individual steps. Finally, a corrected algorithm that takes m steps in the RMVS mapping 
with step size t r and then switches to the WH mapping before taking n steps of size Tw 
can be written as 

£'(mr R + nr w ) = C R l (r R ) o [£ R (r R )] m o C r (t r ) o C^(t w ) o [S w {r w )} n o C w (r w ) (12) 



- 12 - 



where the R subscripts refer to numerical integrations and correctors that perform drifts 
about the Sun (RMVS map) and the W subscripts mark integrations and corrections 
performing drifts about the barycenter (WH map). 

Now that we have an algorithm defined that combines the advantages of both the 
RMVS and WH maps, it is necessary to determine where to switch between maps. As 
one moves further from the planetary region the solar system's potential approaches that 
of a monopole and test particle motion will consequently approach Keplerian motion. At 
. 50 AU, the magnit ude of the quadrupole moment is only 10 -6 of that of the monopole 



(IDuncan et al. 



19871 ) . Through trial and error, we have found the code to behave well if we 
switch maps at a barycentric distance of 300 AU, where the motion of test particles is very 
smooth. (Alternative transition radii are discussed in the next section.) 

In regimes of very smooth test particle motion, extensions to higher order 
approximations will yield large increases in integrator accuracy. As demonstrated in the 
error analysis of the next section, we have found that only a third-order (2-stage) symplectic 
corrector is suitable for our purposes. This is advantageous because the computational 
expense of cor rectors scales d irectly with stage number. The coefficients for this third-order 



corrector are ( lWisdomll2006l ) 



10 



7, h 



1 

72 



7 



1 h - 1 

do = -7, 02 = 7 

5 l, 24 I 



(13) 



where 7 = 10 1 / 2 . 



Simply switching maps and timesteps after a particle crosses the r = 300 AU bo undary 



Quinn et al. 



(119971 ) show 



will produce an integration that is not time-reversible. However, 
that using a time-reversible step-selection routine significantly reduces the error incurred 
from varying timesteps in symplectic integrations. As a result, we use an algorithm that 



attempts to enforce time-reversibility, i.e. the same sequence of step sizes are used if a 



13 



particle is integrated in reverse. To do this, we first drift all particles about the solar system 
barycenter for Tw Only particles that have both initial and final positions further than 300 
AU are integrated using tw and the WH map. All other particles are reset to their initial 
positions and integrated with tr and the RMVS map. It must be mentioned that even 
with this algorithm, time-reversibility is still occasionally violated when the end positions 
of the RMVS, t r integration and the WH, t w drift fall on opposite sides of the r = 300 AU 
boundary. Given the smooth nature of particle motion in this region of the solar system, 



however, this occurrence is re 



over more simplistic routines ( jQuinn et al.lll997| ) 



atively rare, and t his algorithm is still a vast improvement 



Now that we have outlined our numerical method, we look at its performance in the 
next section. 



3. Code Performance 

In this section we look at the performance of our method under a few test cases and 
compare this with other simulation techniques, such as SWIFT's RMVS3 as well as cruder 
multi-stepping methods. We examine numerical errors in energy and angular momentum, 
the effect of passing stars and the Galactic tide on distant orbits, and the change in Jacobi 
constant in the restricted three-body problem. In addition to this, we investigate how 
numerical error changes as we vary T\y and the step/map transition radius, r trans . Last, 
we measure how the computation time of our code scales with particle number. Unless 
specified otherwise, simulations using SCATR are set to tr = 0.548 yrs, tw = 90 = 49.3 
yrs, and r trans = 300 AU. 
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3.1. Numerical Error 

We begin our error analysis with a 1-Gyr integration of a test particle on an orbit 
with semimajor axis a = 300 AU, perihelion distance q = 100 AU, and inclination i = 45° 
under the gravitational influence of the Sun and the four giant planets on their present day 
orbits. Since this particular orbit does not approach the planets but does pass through the 
map/timestep transition radius of our new code twice per orbit, it provides a good measure 
of the numerical error in Keplerian energy and angular momentum. Of course, a particle 
under the influence of the Sun and planets does not actually conserve these quantities, but 
because this orbit stays far from the planets, these quantities change very little and allow 
numerical error to be isolated. 

In Figure |2j we plot the energy error for three different integrators: our new code 
(SCATR), a version of SWIFT RMVS3 with a symplectic corrector, and the normal version 
of SWIFT RMVS3. We see that after 1 billion years the integration using our new code 
has an energy error of less than one part in 10 5 . This error is about 30% larger than that 
of a symplectically corrected integration done entirely in the RMVS map with a constant 
timestep of 0.548 yrs. However, most solar system symplectic integrations do not use a 
corrector, and when we compare the error of SCATR with a normal RMVS3 integration we 
see that although ours shows some secular drift, our error is over an order of magnitude 
smaller than the uncorrected routine. 

An examination of orbital angular momentum evolution for the three integrations yields 
similar results. The change in orbital angular momentum is plotted for each of the three 
simulations in Figure [3J Unlike the energy, angular momentum evolves noticeably during 
our test particle integration. The reason for this is that the higher multipole moments of 
the solar system's potential produce a small secular change in the angular momentum, or 
L, over time. As can be seen in the upper two panels of Figure EJ the angular momentum 
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Fig. 2. — Plots of orbital energy error vs. time for 1-Gyr integrations of a test particle orbit 
with a =300 AU, q =100 AU, and % =45°. Top Panel: This integration is performed with 
our hybrid code, SCATR. Middle Panel: This integration is performed with a version of 
SWIFT RMVS3 outfitted with a third-order symplectic corrector. Bottom Panel: This 
integration is performed with the SWIFT RMVS3 integrator. 
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evolution of the SCATR simulation is nearly identical to the RMVS version that includes a 
corrector. Even though the solar system potential closely resembles a monopole at distances 
beyond ~150 AU, the orbital evolution in our hybrid simulation is much more sensitive to 
the weak effects of the higher order multipoles of the solar system rather than the secular 
effects caused by numerical error, which is encouraging. Lastly, the bottom panel of Figure 
E] demonstrates that the normal version of SWIFT undergoes the same evolution in L, 
but experiences spurious oscillations that are an order of magnitude greater than when a 
corrector is employed. 

It should be noted that with such a large t w (~50 yrs) our sampling frequency of 
the giant planet perturbations is less than some of their orbital frequencies. However, the 
frequency of a perturbation due to a single orbit of a planet is extremely high (as well as 
small in magnitude) compared to the orbital periods of test particles being integrated with 
Tw- Consequently, due to the averaging principle, these weak, high frequency perturbations 
should not affect the long-term evolution of distant test particles. Moreover, the smooth 
change in L and its match to a fixed-step integration in Figure [3] demonstrates that our 
Tw sampling is still capturing the long-term deviations of the solar system potential from 
a monopole (the higher multipole moments) that are critical to the evolution of distant 
orbits. (We tested earlier versions of our code that instead modeled the distant solar 
system potential as just a monopole or a monopole with fixed J2 and J4 terms. However, 
this resulted in the systematic inclusion and exclusion of higher multipole terms at specific 
test particle distances. This produced systematic errors that accumulated quickly for some 
test particle orbital inclinations where these multipole terms significantly impacted orbital 
evolution.) 

Figures H] and [3] illustrate that the numerical error of SCATR is comparable to widely 
used simulation methods over Gyr timescales. However, we have yet to demonstrate how 
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Fig. 3. — Plots of orbital angular momentum vs. time for 1-Gyr integrations of an orbit 
with a =300 AU, q =100 AU, and % =45°. Top Panel: This integration is performed 
with SCATR. Middle Panel: This integration is performed with a version of SWIFT 
RMVS3 outfitted with a third-order symplectic corrector. Bottom Panel: This integration 
is performed with the SWIFT RMVS3 integrator. 
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vital the symplectic corrector is to obtaining the small numerical errors seen in the previous 
two figures. Figure H] illustrates this by comparing the energy error growth of SCATR 
with that of an integration in which a symplectic corrector is not used during timestep 
transitions. In this experiment, we use the two different integrators to integrate the orbits 
of twenty test particles for 1 Gyr. The perihelia and semimajor axes of these orbits are set 
to 100 AU and 250 AU respectively to isolate the numerical energy error and produce the 
highest number of time step transitions per particle (inclinations were chosen randomly). 
In Figure HI we plot the growth of the mean numerical energy error for our suite of particles 
for each of the algorithms. As can be seen in the plot, if the use of a symplectic corrector is 
excluded during timestep transitions, the energy error increases by 2-3 orders of magnitude. 
Thus, it is critical to use a symplectic corrector if we wish to minimize the numerical error 
due to timestep transitions in our algorithm. 

To isolate the numerical error of our integrator from other effects we have purposely 
not included other gravitational perturbations that become important at large distances 
from the Sun, such as passing stars and the galactic tide. Even at orbital distances of a 
few hundred AU, both stellar passages and the galactic tide will drive changes in orbital 
energy and angular momentum that are much larger than those caused by our numerical 
error. When we rerun our previous integrations a nd include a populat ion of passing stars 



(IRickman et al. 



20011 ) . we see in Figure 



20081 ) and a galactic tidal model (ILevison et al. 
[5] that the energy drifts in each simulation are much greater than before, and are caused 
by stochastic impulses from passing stars. Figure |5] also shows that the changes in orbital 
angular momentum are now much larger than previously due to both tidal and stellar 
perturbations. Lastly, the fact that our two integrations experience E and L drifts of similar 
magnitude but finish with different endstates indicates that chaotic evolution under passing 
stars and the galactic tide causes slightly different orbits to diverge on Gyr timescales. 
When these effects are considered, the small errors due to timestep and coordinate changes 
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Fig. 4. — A comparison of the numerical energy error vs. time for an integration using 
SCATR (triangles) and an integration where the use of a symplectic corrector during timestep 
changes is disabled (diamonds) These integrations include the Sun and four giant planets as 
well as 20 test particles. The mean energy error of each suite of particles is plotted. 
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become even less significant. 
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Fig. 5. — Plots of fractional change in orbital energy (dotted line) and angular momentum 
(solid line) vs. time for an orbit with a =300 AU, q =100 AU, and % =45°. These integrations 
include the effects of passing stars and the Milky Way tide. Top Panel: This integration 
is performed with SCATR. Bottom Panel: This integration is done with a version of the 
SWIFT RMVS3 integrator that includes a symplectic corrector. 

Our code was conceived in order to increase the efficiency of Oort Cloud integrations, 
and because forming the Oort Cloud involves repeated scatterings by giant planets, another 
important test is to measure how well our code handles close encounters between test 
particles and massive bodies. To evaluate this, we take the case of the restricted 3-body 
problem. In this simulation, Neptune is placed on a circular orbit around the Sun at 30 
AU and a test particle on an orbit with semimajor axis a = 200 AU, perihelion distance 
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q = 27 AU, and inclination z = 10°. This system is then integrated for 100 Myrs. This 
particular orbit configuration was chosen because it forces the particle to interact strongly 
with Neptune as well as cross the timestep/coordinate transition boundary at r = 300 AU 
twice per orbit. 
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Fig. 6. — Plot of fractional change in Jacobi Constant vs. time for two integrations of a test 
particle under the gravitational influence of a circularly orbiting Neptune mass planet at 30 
AU. The first integration was performed with the SWIFT RMVS3 package (crosses), and 
the second integration was performed with SCATR (diamonds) . 



Two integrations of this system are performed: one using the RMVS3 version of 
SWIFT and the other using SCATR. The results of these integrations are shown in Figure [6] 
where we plot the drift in the Jacobi constant of the test particle with time for each of our 
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integrations. This figure demonstrates that both integration packages experience similar 
levels of numerical error in terms of change in the Jacobi constant. The reason for this 
similar behavior is that the additional error introduced from our timestep and Hamiltonian 
partition change at r = 300 AU is negligible compared to the numerical error introduced 
by close encounters. Thus, our code seems to integrate planetary encounters as well as the 
SWIFT integration package. 

We have now established that the numerical errors of SCATR are smaller than or 
comparable to that of SWIFT's RMVS3 for solar system simulations that use Ty/ = 49.3 
yrs and r trans = 300 AU. However, different degrees of accuracy and computing efficiency 
can be attained if these two parameters are adjusted. For this reason, we now measure 
how the numerical error of our code changes as we vary % and r trans . To do this, we run 
11 different 1-Gyr simulations, each one using a different combination of tw and r trans to 
integrate 20 test particles in the presence of the Sun and giant planets. To set up these 
simulations, we first choose Tw and r tTans . The different values we use for t w are 9.86, 
49.3, 98.6, and 247 yrs, while the chosen values of r trans are 100, 200, and 300 AU. (These 
particular timesteps are chosen because for SCATR to function properly t w must be an 
integer multiple of tr, which we hold fixed at 0.548 yrs.) With tr and r tran s chosen, we 
next assign the same perihelion distance and semimajor axis to all particles in a given 
simulation to model a "worst-case" orbit - one with a short orbital period that crosses our 
chosen r trans but stays far enough from the giant planets to isolate the numerical energy 
error. The perihelion and semimajor axis used for each r tra ns can be found in the legend of 
Figure [7J Finally, we randomly select the remaining orbital elements of the test particles 
(inclination, longitude of ascending node, argument of perihelion, and mean anomaly). 

The results of our 11 simulations are shown in Figured For each simulation we plot 
the mean energy error per particle as a function of r trans and Tw One can see from this 
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Fig. 7. — Plot of mean energy error vs. r w for 11 different 1-Gyr simulations of 20 test parti- 
cles orbits. r trans was varied between 100 (diamonds) , 200 (triangles), and 300 AU (squares). 
The test particle orbits for each r trans value are listed in the figure legend. Simulations 
included the Sun and four giant planets. 
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figure that for a given timestep the numerical error of the integrator grows by over an 
order of magnitude as r trans is moved from 300 AU to 100 AU. This is for two different 
reasons. First, relative to H^epi Hj nt and therefore e are larger at 100 AU than at 300 
AU. Consequently, the accuracy gain of the symplectic corrector is smaller at 100 AU since 
its error scales with e 2 . Second, setting r tra ns = 100 AU enables particles with smaller 
semimajor axes to cross the timestep/map transition than when r trans is set to 300 AU. 
Since these particles have shorter orbital periods, they will accrue numerical error faster 
than their r trans = 300 AU counterparts. 

The second feature one notices in Figure [7] is the general deterioration of integration 
accuracy with increasing tw- Again, there are two reasons for this. As stated previously, 
the error in a symplectic corrector is of order e 2 r 2 , so the error introduced each time Cw is 
calculated will increase with Tjy. In addition, the planetary perturbations to pure Keplerian 
motion are very tiny but still non-zero outside of r trans . While perhaps not important 
compared with the effects of tidal and stellar perturbations at these distances, the planetary 
forces are nevertheless being undersampled by tw- As a result, the error due to Hj nt will 
not be bounded as in a typical symplectic integration. This will produce an additional 
small secular drift in orbital energy. (Note that we did not run a Tw = 247 yrs, r trans = 100 
AU case since t w would be 1/4 of the test particle orbital periods.) 

With the numerical error of our code characterized, we now move on to analyzing its 
efficiency 

3.2. Computing Time 

In this section we measure the reduction in computing time attained by our code. To do 
this we integrate a sample of "generic" Oort Cloud comets, since this is the orbital regime 
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in which we expect our code to perform best. Our chosen comet orbits have semimajor axis 
a between 10 3 AU and 10 5 and a spatial density which scales as a" 7 / 2 . In addition, the 
comets are given isotropic distributions in the other orbital elements. The comets are then 
integrated under the gravitational influence of the Sun and the four giant planets on their 
present day orbits. Each integration is performed for 1 Myr on an Intel Core 2 duo 2.33 
GHz processor twice: once using SCATR (with Tw = 49.3 years and r trans = 300 AU) and 
once with an unaltererd version of the SWIFT RMVS3 package from which our code was 
derived. 




Fig. 8. — Plot of computing time per comet for 1-Myr integrations of Oort Cloud comets 
vs. total number of comets in a simulation. Two sets of integrations are performed: one 
using the SWIFT RMVS3 integrator (squares) and one using our new integrator, SCATR 
(diamonds) . 
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The results of our test integrations are shown in Figure [8] where we plot computing 
time per comet vs. the number of comets in each simulation. As can be seen in this plot, 
our code is significantly faster than RMVS3 for all comet numbers that were tested. In 
fact, for simulations with more than 1000 comets the computing time of our code is more 
than a factor of 30 shorter! The reason that this speed gain depends on the number of 
comets is that in SCATR simulations with small numbers of comets, most of the computing 
time is spent integrating the orbits of the giant planets. Since high semimajor axis orbital 
integrations are so computationally cheap in SCATR, the total computing time does not 
initially change much as the number of comets increases, causing the time per comet to 
fall quickly at first. Only above N > 1000, when the computing load is dominated by 
comet orbit integrations rather than planetary ones, does computing time per comet reach 
a constant value. On the other hand, test particles and planets are always integrated with 
the same time step in RMVS3 simulations, making their computational expense more 
comparable to one another. Consequently, the computing time per comet is essentially 
constant for all of our RMVS3 simulations. 



4. Summary 

The computing time of orbital integrations that transition between regimes of vastly 
different dynamical timescales can be greatly reduced by using different integration step 
sizes and canonical coordinate systems for each regime. However, repeatedly changing 
these two numerical parameters in a symplectic algorithm introduces an accumulation of 
numerical error that can alter the dynamics of the system being modeled. We have shown 
that symplectic correctors can be used to eliminate the bulk of the error introduced during 
timestep and coordinate changes. Implementation of this technique yields an integrator 
that is dramatically faster for simulations of distant test particle orbits. Additionally, it has 
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energy and angular momentum errors comparable to normal symplectic routines, such as 
SWIFT's RMVS3. By decreasing the computing expense of distant orbit integrations by up 
to a factor of ~30, our code has the potential to greatly increase the statistical significance 
of many outer solar system simulations. We will make our code available to any interested 
users, and they should contact us via email to obtain a copy. 
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